﻿using System;
using System.Text;
using MathNet.Numerics.LinearAlgebra.Double;
using MathNet.Numerics.LinearAlgebra.Generic;

namespace innovations.util.exts.mathdotnet
{
    public static class MatrixExt
    {
        public static Vector<double> MeanHorizontally(this Matrix<double> thisMatrix)
        {
            Vector<double> v = new DenseVector(thisMatrix.ColumnCount, 1.0);
            Vector<double> sum = thisMatrix * v;
            return sum * (1.0 / (double)thisMatrix.ColumnCount);
        }

        public static Vector<double> MeanVertically(this Matrix<double> thisMatrix)
        {
            Vector<double> w = new DenseVector(thisMatrix.RowCount, 1.0);
            Vector<double> sum = w * thisMatrix;
            return sum * (1.0 / (double)thisMatrix.RowCount);
        }

        public static Vector<double> SumHorizontally(this Matrix<double> thisMatrix)
        {
            Vector<double> v = new DenseVector(thisMatrix.ColumnCount, 1.0);
            return thisMatrix * v;
        }

        public static Vector<double> SumVertically(this Matrix<double> thisMatrix)
        {
            Vector<double> w = new DenseVector(thisMatrix.RowCount, 1.0);
            return w * thisMatrix;
        }

        public static double Sum(this Matrix<double> thisMatrix)
        {
            var v = new DenseVector(thisMatrix.ColumnCount, 1.0);  // Vector with entries 1.0
            var w = new DenseVector(thisMatrix.RowCount, 1.0);    // Vector with entries 1.0
            return w * thisMatrix * v;
        }

        public static Matrix<double> PointwiseExp(this Matrix<double> m)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount);
            for (int i = 0; i < m.RowCount; i++)
            {
                for (int j = 0; j < m.ColumnCount; j++)
                {
                    temp[i, j] = Math.Exp(m[i, j]);
                }
            }
            return temp;
        }

        public static Matrix<double> PointwisePower(this Matrix<double> m, double power)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount);
            for (int i = 0; i < m.RowCount; i++)
            {
                for (int j = 0; j < m.ColumnCount; j++)
                {
                    temp[i, j] = Math.Pow(m[i, j], power);
                }
            }
            return temp;
        }

        public static Matrix<double> PointwiseSubtraction(this Matrix<double> m, Matrix<double> otherMatrix)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount);
            for (int i = 0; i < m.RowCount; i++)
            {
                for (int j = 0; j < m.ColumnCount; j++)
                {
                    temp[i, j] = m[i, j] - otherMatrix[i, j];
                }
            }
            return temp;
        }

        public static Matrix<double> PointwiseAddition(this Matrix<double> m, Matrix<double> otherMatrix)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount);
            for (int i = 0; i < m.RowCount; i++)
            {
                for (int j = 0; j < m.ColumnCount; j++)
                {
                    temp[i, j] = m[i, j] + otherMatrix[i, j];
                }
            }
            return temp;
        }

        /// <summary>
        /// Replaces one matrix column with another
        /// </summary>
        /// <param name="m">matrix to be changed</param>
        /// <param name="v">replacement vector</param>
        /// <param name="index">index of column in matrix to replace</param>
        /// <returns>updated matrix</returns>
        public static Matrix<double> ReplaceColumn(this Matrix<double> m, Vector<double> v, int index)
        {
            for (int i = 0; i < m.RowCount; i++)
            {
                m[i, index] = v[i];
            }
            return m;
        }

        /// <summary>
        /// Replaces one matrix row with another
        /// </summary>
        /// <param name="m">matrix to be changed</param>
        /// <param name="v">replacement vector</param>
        /// <param name="index">index of row in matrix to replace</param>
        /// <returns>updated matrix</returns>
        public static Matrix<double> ReplaceRow(this Matrix<double> m, Vector<double> v, int index)
        {
            for (int i = 0; i < m.RowCount; i++)
            {
                m[i, index] = v[i];
            }
            return m;
        }

        public static Vector<double> ToVector(this Matrix<double> m, bool vertical = true)
        {
            Vector<double> v = new DenseVector(m.ColumnCount * m.RowCount);
            if (vertical)
                for (int i = 0; i < m.RowCount; i++)
                    if (i == 0)
                        v = new DenseVector(m.Row(i).ToArray());
                    else
                        v = v.Merge(new DenseVector(m.Row(i).ToArray()));
            else
                for (int i = 0; i < m.ColumnCount; i++)
                    if (i == 0)
                        v = new DenseVector(m.Column(i).ToArray());
                    else
                        v = v.Merge(new DenseVector(m.Column(i).ToArray()));
            return v;
        }

        public static Matrix<double> Merge(this Matrix<double> m, Matrix<double> other, MergeMode mergeMode)
        {
            if (mergeMode == MergeMode.Vertical && m.ColumnCount != other.ColumnCount)
                throw new Exception("Matrix dimensions don't match.  Merge can't be completed.");
            if (mergeMode == MergeMode.Horizontal && m.RowCount != other.RowCount)
                throw new Exception("Matrix dimensions don't match.  Merge can't be completed.");
            if (m == null)
                return other;
            if (other == null)
                return m;
            if (mergeMode == MergeMode.Horizontal)
                for (int i = 0; i < other.ColumnCount; i++)
                    m.InsertColumn(m.ColumnCount - 1, other.Column(i));
            else
                for (int i = 0; i < other.RowCount; i++)
                    m.InsertRow(m.RowCount - 1, other.Row(i));
            return m;
        }

        public static string ToTable(this Matrix<double> m)
        {
            StringBuilder stringBuilder = new StringBuilder();
            for (int i = 0; i < m.RowCount; i++)
            {
                for (int j = 0; j < m.ColumnCount; j++)
                {
                    stringBuilder.Append(m[i, j].ToString("#.0000e-000").PadRight(12, ' ') + "\t");
                }
                stringBuilder.Append("\n");
            }
            return stringBuilder.ToString();
        }

        public static Matrix<double> RemoveColumn(this Matrix<double> m, int index)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount - 1);
            for (int i = 0; i < m.RowCount; i++)
                for (int j = 0; j < m.ColumnCount; j++)
                    if (j == index)
                        continue;
                    else if (j < index)
                        temp[i, j] = m[i, j];
                    else
                        temp[i, j - 1] = m[i, j];
            return temp;
        }

        public static Matrix<double> RemoveRow(this Matrix<double> m, int index)
        {
            Matrix<double> temp = new DenseMatrix(m.RowCount, m.ColumnCount - 1);
            for (int i = 0; i < m.ColumnCount; i++)
                for (int j = 0; j < m.RowCount; j++)
                    if (j == index)
                        continue;
                    else if (j < index)
                        temp[j, i] = m[j, i];
                    else
                        temp[j, i - 1] = m[j, i];
            return temp;
        }
    }
}
